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We consider a biased molecular junction subjected to external time-dependent electromagnetic 
field. The field for two typical junction geometries (bowtie antennas and metal nanospheres) is 
calculated within finite-difference time-domain technique. Time-dependent transport and optical 
response of the junctions is calculated within non-equilibrium Green's function approach expressed 
in a form convenient for description of multi-level systems. We present numerical results for a two- 
level (HOMO-LUMO) model, and discuss influence of localized surface plasmon polariton modes on 
transport. 
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I. INTRODUCTION 

Optical properties of structures composed of noble 
metals have long been attracting a considerable atten- 
tion due to unique features of such systems in the 
visible spectrumji^i 3 -^ Recent advances in fabrication 
techniques^ along with a tremendous progress in laser 
technologies opened new venues for application of plas- 
monic materials in biology^ integrated optics nanoscale 
imagining^ 8 - and single molecule manipulation^ Physics 
of surface plasmon phenomenon is relatively simple and 
has long been studied. 10 ' 11 In brief, coherent oscillations 
of conductive electrons in a skin-layer of metal known as 
plasmons are capable of producing strong local electro- 
magnetic (EM) fields in the near-field region. It has been 
reported that such "hot" spots can be localized within 
10 nm or less. This along with a great sensitivity to ini- 
tial conditions and geometry makes plasmonic structures 
so attractable for atom/molecule manipulations. 

A natural combination of nanoplasmonics and molec- 
ular response to the generated field started to appear as 
molecular nanopolaritonicS ) 12 ' 13 which studies molecular 
influence on field propagation, and as a tool for develop- 
ing molecular switches^ The latter utilizes nonadiabatic 
alignment of a molecule on semiconductor surface under 
a tip of scanning tunneling microscope. 

Recent developments in experimental techniques ca- 
pable of measuring optical response of current-carrying 
molecular junction a 15 : 16 lead to theoretical formulations 
suitable for simultaneous description of both transport 
and optical properties of molecular devices. 18 ' 19 

While experimental data are measured in real time, 
theoretical description of both transport and optical re- 
sponse so far has mostly been focused on a steady- 
state description. Time-dependent transport usually is 
treated either within kinetic theor y 20 : 21 or within time- 
dependent density functional approach . 22 i 23 i 24 The for- 
mer generally misses broadening of molecular states due 
to coupling to macroscopic contact a 25 ' 26 : 27 and informa- 



tion on coherence^ 8 - although interesting generalizations 
started to appear^ Limitations of the latter are due to 
absence of developed pseudopotentials and fundamental 
necessity to treat finite (closed) systems (see e.g. Ref. 
for discussion). An alternative approach, based on non- 
equilibrium Green function (NEGF) technique, was ini- 

i — II — ^ii — i 

tially formulated in Refs. I31II32II33I . This approach is a 
natural choice for description of open non-equilibrium 
systems. Moreover it provides possibility to describe re- 
sponse of a molecular junction initially under bias to ex- 
ternal time-dependent perturbation (e.g. laser field). 

Here we consider influence of external field specific for 
particular geometry on transport properties and optical 
response of molecular junction. While formulation of 
time-dependent transport within NEGF is general ; 31 : 32 
all the applications so far were restricted to resonant 
single level models only. We propose a variant of the 
scheme capable of dealing with many-level systems. The 
exact calculations are compared to adiabatic pumping 
regime, frequent in the literature on time-dependent 
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were at the lowest order the problem is re- 



duced to a set of quasi-steady-state solutions with time 
dependent (slow timescale) parameters. Also we gen- 
eralize our previous consideration of steady-state opti- 
cal response of current-carrying junction a 36 ' 37 to a time- 
dependent situation. 

The paper is organized as follows. Section HTl presents 
a model of molecular junction. Section Mil describes 
methodology of EM field calculation. Section IIVI de- 
scribes methodology for simulating transport through 
molecular junction subjected to external time-dependent 
field. Adiabatic pumping version is discussed in section 
IVl Numerical results are presented in section fVTl Section 
IVIII concludes. 



II. MODEL 

We consider a two-level system e± t 2, representing high- 
est occupied (HOMO) and lowest unoccupied (LUMO) 
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molecular orbitals (or ground and excited states in the 
many-body language), coupled to two macroscopic elec- 
trodes L and R. The electrodes are considered to be 
each in its own equilibrium with electrochemical poten- 
tials fiL and [in, respectively. We assume that the driv- 
ing (laser field) frequency is smaller than the plasma fre- 
quency, so that usual division of the junction into non- 
equilibrium molecule coupled to free electron reservoirs 
(metallic contacts) is relevant (for a thorough discussion 
of the assumptions see Ref. l3ll ). Local field at the po- 
sition of the molecule is calculated within finite differ- 
ence time domain technique (see section HTT1 for details), 
and is assumed to be an external time-dependent driv- 
ing force causing (de)excitation in the molecule. Fol- 
lowing Ref. in addition to charge transfer between 
contacts and molecule we introduce also energy transfer 
(coupling of molecular excitations to electron-hole excita- 
tions in the contacts). Molecular excitations are coupled 
to a bath of free photon modes (accepting modes), which 
serve as a measurement device of molecular optical re- 
sponse. Hamiltonian of the system is 

H=Ho + V (1) 

i=l,2 

+ ^2 ekcjck + y^tt; a o| t Oa (2) 

ke{L,R} a 
i=l,2;k£{L,R} 

v= £ (y&clcvQdi + vjfr&ckQd*) (3) 

k^k'£{L,R} 




Here d\ (di) and c k (c&) are creation (annihilation) op- 
erators for an electron in the state i of the molecule and 
state k of the contact, respectively. a) a (a a ) is creation 
(annihilation) operator for a photon in the state a, E(t) 
is external time-dependent field, and jlij =< i\fl\j > is 



matrix element of the molecular (vector) dipole operator 
between states i and j of the molecule (i, j = 1,2). We 
assume flu = fl 22 = (or alternatively one can think 
about these contributions being included into definition 
of the state energies £1.2). V et and V en are matrix ele- 
ments for electron and energy transfer between molecule 
and contacts, and V p represents optical response of the 
molecule. 

Below we consider two approaches to transport and op- 
tical response simulations within the model: exact solu- 
tion of the time-dependent Dyson equation and adiabatic 
pumping regime. The former is similar to the procedure 
described in Refs. [3lll32ll33l however it is presented in a 
form convenient for treating a multi-level molecular sys- 
tem (see section IIVI for discussion) . The latter assumes 
that E(t) can be represented as a product of an oscilla- 
tion of frequency ljq with a slowly varying in time (on 
the timescale of 0J0) envelope F{t). In the spirit of the 
Born-Oppcnhcimer approximation F{t) is considered as 
a parameter when solving electronic part of the prob- 
lem. In this case the form of molecule-field interaction 
becomes (within rotating wave approximation) 

- [md\d 2 e iU0t + /kiddie-*"'*) F(t) (4) 

Details of the approach are presented in section fVl 

As usual, we treat the perturbation V, Eq.©, at the 
second order and within noncrossing approximation.— 
Self-energy due to energy transfer (on the Keldysh con- 
tour) is2I 

S e ™(ri,r 2 )= 2J \Vkk'\ 2 gkiT2,Ti)gk'iTi,T 2 ) 

k^k'e{L,R} 

[£22(71, t 2 ) G 2 i(ti,t 2 )1 , s 

_Gi 2 (ri,r 2 ) Gii(n,T 2 )J [b) 

where Gij are molecular Green functions in the lowest 
order of expansion associated with the Hamiltonian Hq, 
Eq.([2|, and are Green functions of free electrons in 
the contacts. Self-energy due to coupling to photon bath 



£ p (Ti,T 2 )=Era 



i-F 1 a(r 2 ,n)G 22 (Ti,r 2 ) 
Sin^J^dt'F^h-t^it') 



8ir 1 ,r 2 )J t ! oo dt'p 12 it')F^t'-t 1 ) 

iF Q (Tl,T 2 )Gll(Tl,T 2 ) 



(6) 



where F a is Green function for free photon and Pijit) = 
— iG^jit, t) is non-equilibrium reduced density matrix. 

Below we discuss methods for calculating external field 
for different geometries, and present approaches to cal- 
culate time-dependent current and optical response of 
driven molecular junction. 



III. ELECTROMAGNETIC FIELD 
SIMULATIONS 

Among various numerical techniques that allow one 
to predict optical properties of plasmonic systems the 
finite-difference time-domain approach (FDTD) is con- 
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sidcrcd to be the most efficient and yet rclativeiy sim- 
ple. FDTD yields data in perfect agreement with ex- 
perimental measurements and results obtained within 
other techniques^ We simulate optical response of metal 
structures utilizing FDTD approach, in which Maxwell 
equations are discretized in space and time following 
Yee's algorithm^. Dispersion of dielectric constant of 
metal, e(w), is taken in the form of the Drude model 



e(ui) = e r 



iTlo 



(7) 



with numerical parameters describing silver for the wave- 
lengths of interest e r = 8.26, ui p = 1.76 x 10 16 rad/sec, 
r = 3.08 x 10 14 rad/sec. 

For simulations of open systems, one needs to impose 
artificial absorbing boundaries in order to avoid reflection 
of outgoing EM waves back to the simulation domain. 
Among various approaches that address this numerical 
issue, the perfectly matched layers (PML) technique^ 
is considered to be the most adequate. It reduces the 
reflection coefficient of outgoing waves at the simulation 
region boundary to 1CT 8 . Essentially, the PML approach 
surrounds the simulation domain by thin layers of non- 
physical material that efficiently absorbs outgoing waves 
incident at any angle. We implement the most efficient 
and least memory intensive method, convolution per- 
fectly matched layers (CPML) 41 - absorbing boundaries, 
at all six sides of the 3D modeling space. Through ex- 
tensive numerical experimentation, we have empirically 
determined optimal parameters for the CPML bound- 
aries that lead to almost no reflection of the outgoing EM 
waves at all incident angles. Spatial steps, 5x = Sy = Sz, 
along all axes are fixed at 1 nm to assure numerical con- 
vergence and the temporal step is St = Sx/(2c), where c 
is the speed of light in vacuum. 

Numerical integration of Maxwell equations on a grid 
within the FDTD framework was performed at the local 
ASU home-built supercomputer utilizing 120 processors. 
An average execution time for our codes is around 20 
minutes. 

A particular advantage of the FDTD method is its abil- 
ity to obtain the optical response of the structure (as- 
suming linear response) in the desired spectral range in 
a single run42 The system is excited with an ultra-short 
optical pulse constructed from Fourier components span- 
ning the frequency range of interest. Next, Maxwell's 
equations are propagated in time for several hundred 
femtoseconds and the components of the EM field are de- 
tected at the point of interest (for our purposes we con- 
sider the detection point where a molecule is located). 
Fourier transforming the detected EM field on the fly 
yields intensities that can be easily processed into the 
spectral response. Since we also have access to the field 
components, we can evaluate the intensity enhancement 
relative to the incident field. This provides the capability 
for straightforward evaluation of 'coupling efficiency' of 
our plasmonic structures in the spectral range of interest. 



IV. TIME-DEPENDENT TRANSPORT 

We are interested in calculating time-dependent cur- 
rent and optical response of the junction. Expression 
for the current at the interface K (K = L,R) between 
molecule and contact is^ 4 - 



Tr [£<(t,t 1 )G > (t 1 ,t) + G > (t,h)'Ef < (t 1 ,t) (8) 
-£> (t, h) G < (t 1 ,t) - G<(t, h) E K (h,t)] 

where is self-energy due to coupling to contact K 
[£|Kn,r 2 )] . . = V lk g k (T llT2 )V k3 (9) 

and r, a, <, > are retarded, advanced, lesser, and greater 
projections respectively. In the wide band limit, when 
escape rate matrix 



[V K {E)) i:j = 2tt J2 V ik V kj S(E - e k ) 



(10) 



fceif 



is assumed to be energy independent and real part of the 
self-energy ([9|) is disregarded, and when time modulation 
is restricted to molecular subspace only, expression (|8|) 
can be reduced to^I 

M*) =Ip(t) -I° K ut (t) (11) 

Ix(t) = - -r / dEf K (E)lmTi [T K A r (t, E)} (12) 



I° K ut (t)=+ -ReTr [T K p{t)} 



(13) 



where fx{E) is Fermi-Dirac distribution in contact 
K and A r (t,E) is time-dependent (one-sided) Fourier 
transform of the retarded Green function G r (t, t'). 

A r (t,E)= [ di'e <B <*- t ' ) G r (t,t') (14) 

J — oo 

In the absence of time-dependent driving A r (t,E) re- 
duces to usual Fourier transform for retarded Green func- 
tion G r (E) = [S-Ho-S 1 '^)]" 1 . In general S r has con- 
tributions (additive within noncrossing approximation) 
from all the processes involved. p{t) in (| 13|) is reduced 
density matrix 



P (t) = -iG<{t,t) 



(15) 



Lesser and greater Green functions are calculated from 
the time dependent Dyson equation 

G > ' < (t,0 = f dh f dt 2 e- iE ^-^ 
■>>,</ 



X 



A r (t 1 ,E)-£>><(E)A a (t 2 ,E) (16) 
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where 



Al j (t,E) = A r ji (t,E) 



(17) 



and A r (t, E) is denned in Eq. (fT4"|) . 

Contrary to our previous consideratio n 3 ^ 37 optical re- 
sponse of molecular junction is calculated as a true pho- 
ton flux into modes {a}, rather than corresponding elec- 
tronic current between molecular orbitals. We start from 
general expression for time-dependent photon flux into 
mode a (the derivation follows the corresponding proce- 
dure 1 for electronic current, the latter can be found in e.g. 
Ref. Hi 



J a (t) =i < 4 (t)a a (t) >= \V%\ 2 fjt, 

x [f< (t, h)g> (h,t) + g> (t, h )f< (h , t) 
- f> (t, h)g< (h ,t)-g< (t, h)F> (t! , t)] 

Here g is two-particle Green function 



(18) 



Q(r,r') 



(19) 



where D = d[d,2 is molecular de-excitation operator. For 
empty accepting mode a expression (fig)) reduces to 



J a (t) = 



\V p \ 2 



dt 1 e iu '^ tl -^g < (t 1 ,t) (20) 



As in Ref. [37J we approximate the two-particle Green 
function by zero-order (in interaction) expression 

g<(t 1 ,t) » -ih[G^{tM)Giz 2 {ti,t) - pi2{t) P 2i{ti)\ 

(21) 

Note that if envelope change in time is slow (on the 
timcscale of oj a ) second term on the right of (f2Tj) can be 
safely disregarded. In this case expression (|20"|) becomes 
equivalent to approximate expression used in Ref. l37l 
Below we calculate frequency resolved 



J(oj,t) ^^J a (t)S(uj - uj a ) 



(22) 



irh 



7Q HRe / dhe^-^G^t^G^i^t) 



and total 



Jtot(t) 



dujj(u!, t) 



photon fluxes. Here 7 Q (w) = 27r^ Q 5(w 
simulations we use^ 

7 Q (w) = r]U)e~ u/,JJe 



(23) 
w a ),and in 

(24) 



To calculate time-dependent charge, Eq.fTT]), and pho- 
ton, Eq . (|2"0"1) . fluxes one needs time-dependent Fourier 



transform of retarded Green function, Eq. (ri"4|) . The 
Dyson equation for retarded Green function is 

i?- t -H(t))G r (t,t>) (25) 

/+oo 
dt x E r (t - h) G r (t u t') = S(t - t') 
-oo 

Its one-sided Fourier transform leads to equation for 
A r (t, E) in the form 



i^-[H (t)-E] 



dh Y l r (t-t 1 )A r {t ll E) 



(26) 



We consider situation when time-dependent external field 
is applied at time to to a biased molecular junction ini- 
tially at steady-state. In this case differential equation 
(|26p can be solved numerically starting from known ini- 
tial condition A r (t , E) = G Q (E) = [^-Hg-S r (^)]- 1 . 

Alternatively, splitting Ho(t) into time- independent 
H§ and time-dependent Hg(t) parts (average over time 
of the time-dependent part can be included into the time- 
independent Hamiltonian) , one can rewrite Dyson equa- 
tion (|2"5"|) in the integral form 

G r (t, t') = G r Q (t-t')+ f dtx G r (t-tx)H (*i)G r (*i, f) 

(27) 

One-sided Fourier transform of (|27|) leads to integral 
equation for A r (t, E) 



A r (t,E) = G r (E) 



(28) 



+ f G r (t - ti)e iE ( t ~ tl ^~H, l (tx)A r (tx, E) 

Jto 

where lower limit of the integral in the right is set to to 
since H (t < to) = 0. Its solution is 



A r (t,E) =\J eff (t,t Q ;E) G r Q (E) 



(29) 



V eff {t,t ;E) =Tcxp 



[dtx G r (t - tx)e tE ^n (tx) 

J to 

(30) 



Effective evolution operator U e y/ can be obtained by 
variety of methods available in the literature (see e.g. 
Ref. |4y and references therein). One of the simplest 
schemes is cumulant (or Magnus) expansiom 47 i 48 ' 49 

Note that although our consideration is restricted to 
the case when time-dependent driving takes place in the 
molecular subspacc only, generalization to driving in the 
contacts or at the molecule-contact interface is straight- 
forward. 



V. ADIABATIC PUMPING REGIME 

When time evolution of an envelope F(t), Eq.(T3]), is 
slow on the timcscale of the field frequency too- consid- 



cration of the time dependent transport is simplified by 
invoking adiabatic assumption (treating F (t) as a param- 
eter). 

We start with Hamiltonian ([I]) in which interaction 
with driving field is written in the form presented in 
Eq.Q. Transforming the Hamiltonian into rotating 
frame of the field^^i 



8 

w 



6 = — (ni - n 2 ) 



(31) 
(32) 



where n, = (i = 1,2), leads to 



H =H + V (33) 
#o = E - \U\2d\d2 + A*2i4^i) 

i=l,2 

+ E £ fc4 £ * + '^2^ a a [ a a a (34) 

feG{L,fl} « 

+ E (^v* ( - i)W/2 +ff. c .) 

i=l,2;fce{i,R} 

? = E (v^4^4dV Wot + #.c.) (35) 

fc^fc'e{L,.R} 



+ 



E (yS&aQdie**** + H.c. 



where 



(-l)*w /2 



(36) 



Within rotating wave approximation only diagonal el- 
ements of the self-energy due to coupling to the con- 
tacts (electron transfer) £ et , Eq.©. and self-energy due 
to coupling to electron-hole excitations (energy transfer) 
£ e ™, Eq.©, survive 

Sf!(Ti,T2) =S?*(n,r 2 ) e i (- 1 ) < ^( tl - i2 )/ 2 (37) 
S|f(ri,r 2 ) =S|f(n,r 2 )e i (- 1 ) <W0 ( tl - i2 ) (38) 

For self-energy due to coupling to photon bath S p , 
Eq.©, we neglect non-diagonal terms, since they con- 
tribute to retarded (advanced) projection only and cou- 
pling to the bath is assumed to be small relative to cou- 
pling to the contacts. The self-energy becomes diagonal 



S^(ri,T 2 ) = Ef 4 (ri,r 2 )e 



i(— l)*U) (tl— *2) 



(39) 



Resulting Green functions G(ti,t 2 ) depend paramet- 
rically on slow time variable t = (ti + i 2 )/2 through 
time dependence of the envelope F(t), Eq. (IM)) . Trans- 
forming to Wigner coordinates, taking Fourier transform 
in the relative coordinate ti — t 2 , and using gradient 



expansion^ leads to the following expressions for charge 
7 / x v 2 - i n f + °° dE 

M*) = E^ry oo ^ (40) 

5 n S^(£;) d n G>(t,E) d n Tr^{E) d n G<(t,E) 



Tr 



as™ 

and photon 



at™ 



at" 



JaW = l^| 2 E 



■n-f-m 



n.m— 



2™+ m n!m! 



2^r 



(41) 



nn £m \ / am am 



dt n dE r > 



fluxes. Eqs. (j4"0")) and (|4"Tj) are main results of this section. 
They are to be compared with general expressions © and 
(|2U|) . respectively. 



VI. NUMERICAL RESULTS 

We calculate time-dependent transport and optical re- 
sponse by envoking Runge-Kutta scheme with adaptive 
stcpsize control^ 2 - to solve numerically system of differen- 
tial equations ([26]) . 
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FIG. 1: (Color online) Current on the left, II-, and right, la, 
interfaces vs. time for single level model. Numerical results 
(dashed line, red) are compared to analytical expression (solid 
line, blue). Also shown is sum of the currents, II + Ir at the 
two interfaces (dotted line, black). See text for parameters. 

To check accuracy of our numerical approach we start 
from a test calculation for a single level model. Ana- 
lytical solution is available for the latter^ In a biased 
junction (}i L = 1 eV and = —1 eV) the level is set 
below both chemical potentials (eq = —2 eV), so that ini- 
tially the level is occupied and current through the junc- 
tion is negligible (escape rates are = Tr = 0.2). At 
time to position of the level is shifted to eV (stcplike 
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FIG. 2: (Color online) Results of FDTD simulations. Left 
panel shows intensity enhancement as a function of the inci- 
dent wavelength (in nm) in logarithmic scale for two spheres 
of 20 nm in diameter with a gap of 10 nm (solid line, black) 
and bowtie antenna with a gap of 10 nm (dashed line, red), 
and 5 nm (dash-dotted line, blue). Top right inset represents 
steady-state intensity enhancement distribution in logarith- 
mic scale for two spheres system at the resonant wavelength 
of 368.202 nm. Lower right inset shows intensity distribution 
for the bowtie antennas with a gap of 5 nm at 602.647 nm. 



modulation). Here and below we assume Fermi distri- 
butions in the leads corresponding to room temperature 
T = 300 K. Figure [T] presents transient current at the 
two interfaces (direction from contact into the system is 
taken to be positive for both currents) calculated numer- 
ically (dashed line) and with analytical solution (solid 
line). Also shown is sum of the currents at the two in- 
terfaces (dotted line). Outflux of electrons from initially 
fully populated level into the right contact leads to ring- 
ing effect. Eventually the current achieves steady-state. 
Our numerical procedure is seen to give good correspon- 
dence with the analytical result. Below we use similar 
parameters for calculation of time-dependent response of 
the two-level system. 

We consider two geometries of a junction: a bowtie an- 
tenna like electrodes and electrodes in the form of metal- 
lic spheres. Large single-molecule fluorescence measure- 
ments were reported recently for the former^ The latter 
(molecule between two metallic nanoparticles) is custom- 
ary in experimental setups. 

Both structures are excited by a plane wave polarized 
along the axis of symmetry (i.e. along the axis connect- 
ing centers of two spheres, for instance). The electric 
field amplitude is then detected as a function of time. 
Recorded amplitudes are Fourier transformed and nor- 
malized with respect to the incident field amplitude lead- 
ing to enhancement as a function in the frequency do- 
main. 
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FIG. 3: (Color online) Comparison of exact numerical solu- 
tion (solid line, red) to adiabatic approximation (dashed line, 
blue) for the two-level (HOMO-LUMO) model. Shown are 
(a) levels populations and (b) current at the left interface vs. 
time. See text for parameters. 



Results of our simulations for both geometries are pre- 
sented in Fig. [2] showing intensity enhancements in the 
main panel. As expected bowtie structures result in no- 
ticeably higher enhancements reaching 630 centered at 
A = 600 nm for a bowtie antenna with a gap of 5 nm. 
Two spheres also show significant enhancement of 55 
around A = 370 nm. We note that the bowtie antenna 
in comparison to two spheres system exhibits two reso- 
nances. The "blue" resonance located at low wavelength 
corresponds to rod lightning effect with high enhance- 
ment localized primarily at the edges of each triangle. 
This feature disappears from the spectrum once sharp 
corners are replaced with smooth edges42 Top and bot- 
tom insets show intensity enhancement distributions at 
resonant conditions for the two spheres and bowtie an- 
tennas, respectively. We place molecular junction in the 
hot spot regions. 

Figure [3^ shows time-dependent populations of molec- 
ular junction driven by external electromagnetic field for 
the ground, n\, and excited, n 2 , states. Time-dependent 
current at the left interface, 1^. is shown in Fig. [3Jd. Pa- 
rameters of the calculation are T = 300 K, e\ = — 1 eV, 
e 2 = 1 cV, [T K ] mm = 0.1 cV and [T K } 12 = [T K } 21 = 
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FIG. 4: (Color online) The two-level (HOMO-LUMO) model. 
Shown are (a) current and (b) total optical response, (|23[) . 
vs. time for bowtie nanoantennas (the strongest signal, blue) 
and two spheres junction geometries. In the latter case the 
response is calculated for two positions of the molecule in 
the junction: in the middle between the spheres (the weakest 
signal, red) and closer to one of the spheres (intermediate sig- 
nal, (a) white silhouette and (b) solid line, black). Figure (c) 
shows contour map of optical flux, (|22)) . for bowtie geometry 
vs. outgoing frequency and time. See text for parameters. 



(m = 1,2 and K = L,R). For interaction with elec- 
tromagnetic field we take fiEo = 0.005 eV, where Eq is 
amplitude of the external laser field before enhancement. 
Bias V is applied symmetrically ^l,b, = Ep ± eV/2, and 
the Fermi energy is Ep = 0. Results presented in Fig. [3] 



are obtained for bowtie geometry with 10 nm gap at bias 
V = 2 V. Exact numerical calculation (solid line) is com- 
pared with adiabatic approximation data (dashed line). 
One sees, that the adiabatic approximation for realis- 
tic parameters provides qualitatively correct results. It 
misses however delay (memory) effects and overestimates 
response signal. Electromagnetic pulse depletes ground 
state and populates excited state, which for the chosen 
bias leads to increase of current through the junction due 
to increase in transmission of the excited state channel 
(see also Fig. [5] below). 

We compare response of the two molecular junction 
geometries in Fig. 3J Bowtie geometry provides stronger 
local enhancement, and consequently stronger molecular 
response. In the case of spherical nanoparticles we con- 
sider two possible positions of molecule between the elec- 
trodes: symmetric and asymmetric (3 nm shift from the 
center, where the field enhancement for the geometry is 
strongest). These yield weakest and intermediate signal, 
respectively. Note, that it is natural to expect that local 
field enhancement is stronger for a structure with uneven 
surface. Fig. [4£i presents time-dependent current for the 
three cases. Total optical response, Eq. ([2"3"|) . is shown in 
Fig. [3)3. We choose r\ = 5 • 10~ 5 and oj c = 2 eV, other 
parameters are as in Fig. [3] Note much more sensitive 
character of optical response to resonant conditions. It 
results from our choice of 7 a (w), Ea. p4]l . so that most 
of the electronic excitation contributes to current. While 
the choice is arbitrary, it indicates importance of the en- 
vironment (bath spectral density). Fig. [4J; shows time- 
dependent optical spectrum, Eq. ([22|) . for the bowtie ge- 
ometry. The signal follows (with a delay) the pulse of 
the external field. Asymmetric character of the spec- 
trum relative to resonance, uj a = 2 eV, stems from over- 
lap of Lorentzians (levels boradening due to coupling to 
the contacts) centered on ground and excited states. 




50 100 150 

t(fs) 

FIG. 5: (Color online) Current vs. time for the two-level 
(HOMO-LUMO) model calculated at two different biases. See 
text for parameters. 

Figure [5] shows time-dependent current response to ex- 
ternal driving at two different constant biases. The cal- 
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culatipn is done for bowtie geometry with a gap of 10 nm, 
param&ers-Wc the same as in Fig. [3] For pre-resonant 
bias, V = 1.8 V, optical excitation is effective in deplet- 
ing the ground and populating the excited states of the 
mojDi^i^ijQQrich results in increased current through both 
channels. At post-resonant bias, V = 2.2 V, the charge 
transfer channels are open. Here optical excitation con- 
tributes mostly to decrease in conductance of the ground 
staxc'^ncrappcarancc of leakage current to the left con- 
tact in the excited state. This leads to overall decrease in 
current through the junction (see also discussion below). 

2.8e-06 




50 100 150 

t(fs) 




FIG. 6: (Color online) Role of energy transfer process. Shown 
are (a) total optical response vs. time with (dotted line, red) 
and without (solid line, blue) electron-hole excitations and 
(b) difference between current calculated with and without 
electron-hole excitations vs. time and bias. Calculations are 
performed within adiabatic approximation scheme. See text 
for parameters. 

Calculations so far disregarded influence of both en- 
ergy transfer, Eq.©, and external photon bath, Eq.©, 
(except its contribution to optical rate) on electronic dis- 
tribution in the molecule. While the latter can indeed be 
disregarded due to smallness of the reasonable coupling 



parameter (see Ref. for discussion), the former can 
make a difference. Here we illustrate influence of energy 
transfer process on time-dependent response of the junc- 
tion within adiabatic approximation (full numeric calcu- 
lation is straightforward but time-consuming). Figure [Bp, 
shows total optical response calculated with (dashed line) 
and without (solid line) energy transfer included. Calcu- 
lation is done for bowtie geometry with 10 nm gap at 
pre-resonant constant bias V = 1.8 V. Other parameters 
are as in Fig. [3] As expected, energy transfer diminishes 
optical response of the junction, since both energy trans- 
fer from molecule to contacts and fluorescence compete 
for the same excess electronic population in the excited 
state. Current change upon including electron-hole exci- 
tations into consideration is more interesting. Interplay 
between channel blocking and resonant pathways for elec- 
tron transfer may lead to increase in current through the 
junction as is illustrated in Fig.[Bb). This effect is similar 
to the situation presented in Fig. [Sj 



VII. CONCLUSION 

We consider a two-level (HOMO-LUMO) model of 
molecular junction driven by external time-dependent 
laser field. Finite difference time domain technique is 
used to calculate field distribution for two junction ge- 
ometries. Resulting local field at the molecule is con- 
sidered to be the driving force. We assume that the 
junction is initially in a nonequilibrium steady-state re- 
sulting from applied constant bias. At time to driving 
force (laser pulse) starts to influence the system. Time- 
dependent transport (charge flux through the junction) 
and optical response (photon flux from the molecule into 
accepting modes) are calculated for a set of geometries 
and applied biases. We rewrite a nonequilibrium Green 
function technique for time-dependent calculation in a 
form convenient for treating many-level molecular sys- 
tems. Results of the simulations within the approach 
are compared to approximate scheme for an adiabatic 
pumping regime. Note that while our present considera- 
tion is restricted to driving force applied to the molecule 
only, generalization of the approach to situations of time- 
dependent bias and/or coupling between molecule and 
contacts is straightforward. Extension of the consider- 
ation to realistic molecular devices, taking into account 
time-dependent non-equilibrium distribution in the con- 
tacts and spatial profile of the field, and considering in- 
terplay of time-dependencies of bias and laser field are 
goals of future research. 
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